home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
AOL File Library: 2,801 to 2,900
/
aol-file-protocol-4400-2801-to-2900.zip
/
AOLDLs
/
C++ Files Library
/
Graphic Gems I, II & III (C_C++)
/
Graphics Gems C Code.sea
/
GemsIII
/
partition3d
/
partition.c
< prev
next >
Wrap
Text File
|
1992-06-16
|
12KB
|
373 lines
/* partition.c: module to partition 3D convex face with an arbitrary plane.
* Copyright (c) Norman Chin
*/
#include <stdio.h> /* NULL */
#include <assert.h> /* assert() */
#include <malloc.h> /* malloc() */
#include <math.h> /* fabs() */
#include "GraphicsGems.h"
#include "partition.h"
/* local functions */
static VERTEX *findNextIntersection(/* VERTEX *vstart, double aa,
double bb, double cc, double dd,
double *ixx, double *iyy, double *izz,
SIGN *sign */);
static SIGN anyEdgeIntersectWithPlane(/* double x1, double y1, double z1,
double x2, double y2, double z2,
double aa, double bb, double cc,
double dd,
double *ixx, double *iyy, double *izz
*/);
static FACE *createOtherFace(/* FACE *face, VERTEX *v1, double ixx1,
double iyy1, double izz1, VERTEX *v2,
double ixx2, double iyy2, double izz2 */);
static SIGN whichSideIsFaceWRTplane(/* FACE *face, double aa, double bb,
double cc, double dd */);
static VERTEX *allocVertex(/* double xx, double yy, double zz */);
static FACE *allocFace(/* FACE *face, VERTEX *vhead */);
/* Partitions a 3D convex polygon (face) with an arbitrary plane into its
* negative and positive fragments, if any, w.r.t. the partitioning plane.
* Note that inputFace is unusable afterwards since its vertex list has been
* parceled out to the other faces. It's set to null to avoid dangling
* pointer problem.
*/
void partitionFaceWithPlane(aa,bb,cc,dd,inputFace,faceOn,faceNeg,facePos)
double aa,bb,cc,dd; /* partitioning plane's coefficients */
FACE **inputFace; /* face to be partitioned */
FACE **faceOn; /* returns face embedded in plane, if any */
FACE **faceNeg; /* returns face on negative side, if any */
FACE **facePos; /* returns face on positive side, if any */
{
VERTEX *v1, *v2;
double ixx1,iyy1,izz1, ixx2,iyy2,izz2;
SIGN signV1, signV2;
*faceOn= *faceNeg= *facePos= NULL_FACE;
/* find first intersection */
v1= findNextIntersection((*inputFace)->vhead,aa,bb,cc,dd,
&ixx1,&iyy1,&izz1,&signV1);
if (v1 != NULL_VERTEX) {
assert(signV1 != ZERO);
/* first one found, find the 2nd one, if any */
v2= findNextIntersection(v1->vnext,aa,bb,cc,dd,
&ixx2,&iyy2,&izz2,&signV2);
/* Due to numerical instabilities, both intersection points may
* have the same sign such as in the case when splitting very close
* to a vertex. This should not count as a split.
*/
if (v2 != NULL_VERTEX && signV1 == signV2) v2= NULL_VERTEX;
}
else v2= NULL_VERTEX; /* No first intersection found,
* therefore no second intersection.
*/
/* an intersection? */
if (v1 != NULL_VERTEX && v2 != NULL_VERTEX) { /* yes, intersection */
FACE *newOtherFace;
/* inputFace's vertex list will be modified */
newOtherFace= createOtherFace(*inputFace,v1,ixx1,iyy1,izz1,
v2,ixx2,iyy2,izz2);
/* return split faces on appropriate lists */
if (signV1 == NEGATIVE) {
*faceNeg= *inputFace;
*facePos= newOtherFace;
}
else {
assert(signV1 == POSITIVE);
*faceNeg= newOtherFace;
*facePos= *inputFace;
}
}
else { /* no intersection */
SIGN side;
/* Face is embedded or wholly to one side of partitioning plane. */
side= whichSideIsFaceWRTplane(*inputFace,aa,bb,cc,dd);
if (side == NEGATIVE)
*faceNeg= *inputFace;
else if (side == POSITIVE)
*facePos= *inputFace;
else {
assert(side == ZERO);
*faceOn= *inputFace;
}
}
/* inputFace's vertex list has been parceled out to other lists so
* set this to null.
*/
*inputFace= NULL_FACE;
/* If a face is embedded in plane then there must not be any faces on
* either side.
* Otherwise the face isn't embedded in plane so there must be at least
* one face on either side or both.
*/
assert( (*faceOn != NULL_FACE) ?
(*faceNeg == NULL_FACE &&
*facePos == NULL_FACE) :
(*faceNeg != NULL_FACE ||
*facePos != NULL_FACE) );
} /* partitionFaceWithPlane() */
/* Finds next intersection on or after vstart.
*
* If an intersection is found,
* a pointer to first vertex of the edge is returned,
* the intersection point (ixx,iyy,izz) and its sign is updated.
* Otherwise a null pointer is returned.
*/
static VERTEX *findNextIntersection(vstart,aa,bb,cc,dd,ixx,iyy,izz,sign)
VERTEX *vstart; /* where to start searching on vertex list */
double aa,bb,cc,dd; /* plane's coefficients */
double *ixx,*iyy,*izz; /* returned intersection point */
SIGN *sign; /* returned classification of intersection point */
{
VERTEX *vtrav;
/* for all edges starting from vstart ... */
for (vtrav= vstart; vtrav->vnext != NULL_VERTEX; vtrav= vtrav->vnext) {
if ((*sign= anyEdgeIntersectWithPlane(vtrav->xx,vtrav->yy,vtrav->zz,
vtrav->vnext->xx,vtrav->vnext->yy,
vtrav->vnext->zz,aa,bb,cc,dd,
ixx,iyy,izz))) {
return(vtrav);
}
}
return(NULL_VERTEX);
} /* findNextIntersection() */
/* Memory allocated for split face's vertices and pointers tediously updated.
*/
static FACE *createOtherFace(face,v1,ixx1,iyy1,izz1,v2,ixx2,iyy2,izz2)
FACE *face; /* face to be split */
VERTEX *v1; /* 1st vertex of edge of where 1st intersection was found */
double ixx1,iyy1,izz1; /* 1st intersection */
VERTEX *v2; /* 1st vertex of edge of where 2nd intersection was found */
double ixx2,iyy2,izz2; /* 2nd intersection */
{
VERTEX *i1p1, *i2p1; /* new vertices for original face */
VERTEX *i1p2, *i2p2; /* new vertices for new face */
VERTEX *p2end; /* new vertex for end of new face */
VERTEX *vtemp; /* place holders */
register VERTEX *beforeV2; /* place holders */
FACE *newFace; /* newly allocated face */
/* new intersection vertices */
i1p1= allocVertex(ixx1,iyy1,izz1);
i2p1= allocVertex(ixx2,iyy2,izz2);
i1p2= allocVertex(ixx1,iyy1,izz1);
i2p2= allocVertex(ixx2,iyy2,izz2);
/* duplicate 1st vertex of 2nd list to close it up */
p2end= allocVertex(v2->xx,v2->yy,v2->zz);
vtemp= v1->vnext;
/* merge both intersection vertices i1p1 & i2p1 into 1st list */
if (ISVERTEX_EQ(i2p1,v2->vnext)) { /* intersection vertex coincident? */
FREEVERTEX(i2p1); /* yes, so free it */
i1p1->vnext= v2->vnext;
}
else {
i2p1->vnext= v2->vnext; /* attach intersection list onto 1st list */
i1p1->vnext= i2p1; /* attach both intersection vertices */
}
v1->vnext= i1p1; /* attach front of 1st list to intersection vertices */
/* merge intersection vertices i1p2, i2p2 & p2end into second list */
i2p2->vnext= i1p2; /* attach both intersection vertices */
v2->vnext= i2p2; /* attach 2nd list to interection vertices */
if (vtemp == v2) {
i1p2->vnext= p2end; /* close up 2nd list */
}
else {
if (ISVERTEX_EQ(i1p2,vtemp)) { /* intersection vertex coincident? */
FREEVERTEX(i1p2); /* yes, so free it */
i2p2->vnext= vtemp; /* attach intersection vertex to 2nd list */
}
else {
i1p2->vnext= vtemp; /* attach intersection list to 2nd list */
}
/* find previous vertex to v2 */
for (beforeV2= vtemp; beforeV2->vnext != v2; beforeV2= beforeV2->vnext)
; /* lone semi-colon */
beforeV2->vnext= p2end; /* and attach it to complete the 2nd list */
}
/* copy original face info but with new vertex list */
newFace= allocFace(face,v2);
return(newFace);
} /* createOtherFace() */
/* Determines which side a face is with respect to a plane with coefficient
* (aa,bb,cc,dd).
*
* However, due to numerical problems, when a face is very close to the plane,
* some vertices may be misclassified.
* There are several solutions, two of which are mentioned here:
* 1- classify the one vertex furthest away from the plane, (note that
* one need not compute the actual distance) and use that side.
* 2- count how many vertices lie on either side and pick the side
* with the maximum. (this is the one implemented).
*/
static SIGN whichSideIsFaceWRTplane(face,aa,bb,cc,dd)
FACE *face;
double aa,bb,cc,dd;
{
register VERTEX *vtrav;
double value;
boolean isNeg, isPos;
isNeg= isPos= FALSE;
for (vtrav= face->vhead; vtrav->vnext != NULL_VERTEX; vtrav= vtrav->vnext){
value= (aa*vtrav->xx) + (bb*vtrav->yy) + (cc*vtrav->zz) + dd;
if (value < -TOLER)
isNeg= TRUE;
else if (value > TOLER)
isPos= TRUE;
else assert(ISDOUBLE_EQ(value,0.0));
} /* for vtrav */
/* in the very rare case that some vertices slipped thru to other side of
* plane due to round-off errors, execute the above again but count the
* vertices on each side instead and pick the maximum.
*/
if (isNeg && isPos) { /* yes so handle this numerical problem */
int countNeg, countPos;
/* count how many vertices are on either side */
countNeg= countPos= 0;
for (vtrav= face->vhead; vtrav->vnext != NULL_VERTEX;
vtrav= vtrav->vnext) {
value= (aa*vtrav->xx) + (bb*vtrav->yy) + (cc*vtrav->zz) + dd;
if (value < -TOLER)
countNeg++;
else if (value > TOLER)
countPos++;
else assert(ISDOUBLE_EQ(value,0.0));
} /* for */
/* return the side corresponding to the maximum */
if (countNeg > countPos) return(NEGATIVE);
else if (countPos > countNeg) return(POSITIVE);
else return(ZERO);
}
else { /* this is the usual case */
if (isNeg) return(NEGATIVE);
else if (isPos) return(POSITIVE);
else return(ZERO);
}
} /* whichSideIsFaceWRTplane() */
/* Determines if an edge bounded by (x1,y1,z1)->(x2,y2,z2) intersects
* the plane with the coefficients (aa,bb,cc,dd).
*
* If there's an intersection,
* the sign of (x1,y1,z1), NEGATIVE or POSITIVE, w.r.t. the plane is
* returned with the intersection (ixx,iyy,izz) updated.
* Otherwise ZERO is returned.
*/
static SIGN anyEdgeIntersectWithPlane(x1,y1,z1,x2,y2,z2,aa,bb,cc,dd,
ixx,iyy,izz)
double x1,y1,z1, x2,y2,z2; /* both vertices of edge */
double aa,bb,cc,dd; /* partitioning plane's coefficients */
double *ixx,*iyy,*izz; /* returned intersection point, if any */
{
double temp1, temp2;
int sign1, sign2; /* must be int since gonna do a bitwise ^ */
/* get signs */
temp1= (aa*x1) + (bb*y1) + (cc*z1) + dd;
if (temp1 < -TOLER)
sign1= -1;
else if (temp1 > TOLER)
sign1= 1;
else {
/* edges beginning with a 0 sign are not considered valid intersections
* case 1 & 6 & 7, see text
*/
assert(ISDOUBLE_EQ(temp1,0.0));
return(ZERO);
}
temp2= (aa*x2) + (bb*y2) + (cc*z2) + dd;
if (temp2 < -TOLER)
sign2= -1;
else if (temp2 > TOLER)
sign2= 1;
else { /* case 8 & 9, see text */
assert(ISDOUBLE_EQ(temp2,0.0));
*ixx= x2;
*iyy= y2;
*izz= z2;
return( (sign1 == -1) ? NEGATIVE : POSITIVE);
}
/* signs different?
* recall: -1^1 == 1^-1 ==> 1 case 4 & 5, see text
* -1^-1 == 1^1 ==> 0 case 2 & 3, see text
*/
if (sign1 ^ sign2) {
double dx,dy,dz;
double denom, tt;
/* compute intersection point */
dx= x2-x1;
dy= y2-y1;
dz= z2-z1;
denom= (aa*dx) + (bb*dy) + (cc*dz);
tt= - ((aa*x1) + (bb*y1) + (cc*z1) + dd) / denom;
*ixx= x1 + (tt * dx);
*iyy= y1 + (tt * dy);
*izz= z1 + (tt * dz);
assert(sign1 == 1 || sign1 == -1);
return( (sign1 == -1) ? NEGATIVE : POSITIVE );
}
else return(ZERO);
} /* anyEdgeIntersectWithPlane() */
/* allocates vertex with coordinates (xx,yy,zz) */
static VERTEX *allocVertex(xx,yy,zz)
double xx,yy,zz;
{
VERTEX *newVertex;
newVertex= NEWTYPE(VERTEX); assert(newVertex != NULL_VERTEX);
newVertex->xx= xx; newVertex->yy= yy; newVertex->zz= zz;
newVertex->vnext= NULL_VERTEX;
return(newVertex);
} /* allocVertex() */
/* allocates face with list of vertices, vhead */
static FACE *allocFace(face,vhead)
FACE *face;
VERTEX *vhead;
{
FACE *newFace;
newFace= NEWTYPE(FACE); assert(newFace != NULL_FACE);
*newFace= *face;
newFace->vhead= vhead;
return(newFace);
} /* allocFace() */